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As a prelude to what might be expected as forthcoming breakthroughs in 
finding new approaches toward solving three-dimensional lattice models in 
the twenty-first century, we review the exact solutions of two lattice models in 
three dimensions obtained using the conventional combinatorial and transfer 
matrix approaches. 
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1 Introduction 



The twenties century has witnessed the birth and growth of the field of ex- 
actly solvable models in statistical mechanics. It began with the Lenz-Ising 
model of Ising[|J and its solution in two dimensions by Onsager.0 Begin- 
ning in the 1960's the field of exactly solvable models grew phenomenally 
to encompass vertex-type models|3|, [|] and related problems. But most of 
the new soluble models are limited to two dimensions. To be sure, several 
three-dimensional systems have been studied with various degrees of suc- 
cess in recent years. They include the Zamolodchikov model solved by 
Baxter, |7]] its more recent extension by Bazhanov and Baxter, j§| and two 
models solved by this author and collaborators. || pi], [□]] But these solutions 
invariably suffer defects of one kind or another, and the prevailing consensus 
is that the solution of genuine three-dimensional systems must await for ma- 
jor breakthroughs. Before the arrival of the new breakthroughs, however, it is 
instructive to pause, and examine the usefulness of conventional approaches 
when applied to three-dimensional systems. 

There are generally two conventional approaches, namely, a combinato- 
rial approach using graphical methods including Pfaffians and an algebraic 
approach using the method of the transfer matrix, both of which have been 
applied with some success in treating three-dimensional systems. Here, we 
review two such solutions. |9|, |10| It is hoped that, with an understanding 
of what has been done and what cannot be done with the conventional ap- 
proaches, one can then proceed ahead to explore new lines of inquiry toward 
solving genuine three-dimensional lattice models in the twenty-first century. 



2 The combinatorial solution of an 0(— 1) model 
2.1 Definition of the model 

Consider a simple cubic lattice C of M = N\ x N2 x ^3 sites. The partition 
function of an 0(n) model ]I2], on £ is the graph generating function 

Z(n)= £ nV, (1) 

closed polygons 

where the summation is taken over all closed non-intersecting polygonal con- 
figurations on C, I is the number of polygons, and b is the number of edges 
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contained each configuration. We shall restrict to polygons which are di- 
rected in a generally negative to positive direction along the three principal 
axes, hence a special three-dimensional system. Assuming periodic bound- 
ary conditions, the lines will then loop around the lattice, with each loop 
cutting one or more times a plane perpendicular to the (1, 1, 1) direction. 
While the evaluation of (H) for general n is still open, we consider here the 
n — — 1 problem and are interested in the per-site partition function, or "free 
energy" , 

/= lim rinZf-l). (2) 



2.2 Equivalence with a dimer problem 

Our strategy is to establish a one-one correspondence between the configu- 
rations in (|ip and dimer configurations on a related lattice £', and that the 
latter can be evaluated as a Pfafnan. To construct £', we split each lattice 
point of £ into two, with one point attached to the three positive axes and 
the other to the three negative axes as shown in Fig. 1. The dimer lattice 
£' is then obtained by reconnecting the split pairs with new inserted edges. 
We shall associate dimer weights z to edges originally in £, and weights 1 to 
the newly inserted edges as shown. 

The one-one correspondence between vertex configurations on £ and 
dimer configurations on £' can be seen by superimposing any given dimer 
configuration d with a standard one, Co, in which all inserted edges are cov- 
ered by dimers. The superposition of the two dimer configurations produces 



a graph of transition cycles, |14[ or polygons, belonging to one of two kinds: 
i) double dimers placed on new edges, and ii) line configurations looping 
around the lattice. In the latter case the line configurations form precisely 
the same polygonal configurations as in ([[]). Conversely, to each polygo- 
nal configuration in ([!]), there exists a unique dimer configuration C, which, 
when superimposed with Co, produces the line configuration in question. 
This establishes the desired correspondence. 

Now there are M sites on £, and this leads us to consider a Pfafnan Pf A 
of dimension 2J\f. The convention of fixing signs to lattice edges is to direct 
edges of £' such that the element A a/3 is positive (negative) if an edge is 
directed from site a to (3 (/3 to a). We choose to direct all edges of £' along 
the positive directions as shown in Fig. 1. Consider a typical term in the 
Pfafnan corresponding to a dimer configuration Cj. The sign of this term 
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relative to the term corresponding to Co is the product of the signs of the 
transition cycles produced by the superimposition of Cj and Co- The rule is 
that each transition cycle carries a sign (— l) m+1 , where m is the number of 



arrows pointing in a given direction when the transition cycle is traversed. [|14]] 
For transition cycles consisting of double dimers we have m = 1 so that the 
sign is always positive. For other transition cycles, or loops, we have always 
m = even, so that each loop carries a factor —1. As a consequence, the 
relative sign of the term is (— 1)^ where £ is the number of loops. This 
establishes that, with our choice of signs for A a p, the Pfaffian A is precisely 
the dimer generating function (2) with n = — 1. 



2.3 The free energy 

The Pfaffian Pf A is the square root of a 2J\f x 2J\f antisymmetric matrix A. 
Explicitly, we have 

Z(-l) = Pf A = VdetA (3) 

where 



A 







-1 



(g) I N1 <g> I N2 <g> I Na 



z 
-z 



H Nl ® I N2 ® I N3 



+1 



H 



N 2 ^ In 3 + Ini 



In, ® ^ 



(4) 



Here, Ja? is the N x N identity matrix, Hn is an N x N matrix given by 



H 



N 



( 10. 

1. 

0. 

V-i o o . 



o\ 



1 



(5) 



Now, ifjv can be diagonalized with eigenvalues e^, where = 2irn/N, n = 
1, 2, • • • , N. Replacing the H matrices by diagonal ones with its eigenvalues 
in the diagonal, one obtains 



JVi N 2 N 3 r 

z(-i) = n n n 

ni=l «i=l n3=l 



1 + Ei ze 2mn ^ N * 

-1 - Ei-ze" 2 ™ 1 /^ 



1/2 



4 



Ni N 2 N s 

nnn 

ni=l «2=1 U3=l 



i + E 



ze 



2irini/Ni 



i=l 



This leads to the thermodynamic limit of the free energy 
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] p27T f-l'K f-l'K 



i=l 



(6) 



2.4 The critical behavior 



The analysis of the critical behavior of the free energy ([]) is facilitated by 
the use of the integration formula 



1 r 2n 

— / d8ln\A + Be i0 \ 



lnmax{|A|, \B\} 



(7) 



which holds for any complex A and B. We realize bond weights by intro- 
ducing bond energy e, and write without the loss of generality z = e _/3<E < 1, 
where (3 = 1/kT. Then, multiplying the quantity inside the absolute signs 
in (0) by e~ ldl and carrying out the integration over 9%, we obtain 



./ 



0. 



T <T r 



1 



(2tt) 2 

where T c is defined by 



d9o ■ ■ ■ I d6 3 In 

n 



3z 



3=2 



T > T r 



(8) 



(9) 



and the integrations are taken over the regime K specified by z\ 1 +e i(>2 +e i6s | > 
1. Clearly, the free energy is non-analytic in T at T c with the system frozen 
below T c . 

To find out the nature of the transition, we evaluate the energy 



U 



d 



(2tt) 2 V dz 



/ d9\ / d02 lnmaxjx, \x\} 
Jo Jo 



d9 1 / d0 2 



n 



(10) 
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where 

x= l + z(e Wl +e id2 ) } (11) 



and 1Z is the regime z > \x\. In writing out the second line in (|10D, we have 
carried out the integration over # 3 . The last line in (|TUp then follows from 
the fact that x is independent of z so that the regime \x\ > z gives rise to no 
contribution after taking the derivative. 

It is clear that the volume of TZ is small near T c . A detailed analysis [|T0l 
shows that 

U~\T-T C \ (12) 
and hence the critical exponent a = 0. Explicitly, one finds 

U = 0, T<T C 

2e f 2z du _ x ( l-u 2 -z 2 \ 

cos , T >T C (13) 



TC 2 Jz-1 \/Az 2 — u 2 V 2uz 

from which one obtains 

dU 3v3,, ,o, m m , 

° v = ~dT = ^T (ln3) K T^T C + . (14) 

Thus, the specific heat has a cusp singularity at T c +. 

The analysis of this section can be extended to the 0(— 1) model defined 
by (HD in d dimensions in which continuous lines run in a preferred direction 
along all d principal axes of a rf-dimensional hypercubic lattice. In this case 
one again finds the system frozen below T c , but with T c defined by dz = 1. 
One also finds the critical exponent a — (d— 3)/2. 



3 The transfer matrix solution of an interact- 
ing layered dimer system 

3.1 Definition of the model 

Consider a three-dimensional model consisting of K layers of honeycomb 
dimer lattices. The dimers, which carry weights u,v,w along the three hon- 
eycombO edge directions, are close-packed within each layer and, in addition, 
interact between layers. The situation is shown in Fig. 2. For two dimers 
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Table 1: Dimer interacting energy between two dimers incident at the same 
site of adjacent layers. 



layer k — > k + 1 


w 


f 




w 





-2/i/3 


2/i/3 


V 


2/i/3 





-2/i/3 


u 


-2/i/3 


2/i/3 






incident at the same site in adjacent layers, the interaction energy is given 
in Table 1. 

The two-dimensional honeycomb dimer system can be formulated as a 



five- vertex model, namely, an ice-rule model with the weights) 15 



{Ul, Co> 2 , UJ 3 , Co> 4 , LU 5 , Co> 6 } = {0, W, V, It, \/uV, \/uv} . (15) 

The 5-vertex model is defined on a square lattice of size M x N mapping to 
an honeycomb lattice of 2MN sites. |J5] The mapping is such that the edge 
state a = +1 (P = +1) corresponds to the presence of a v (u) dimer. 

The K layers of M x N square lattice form a a simple-cubic lattice C of 
size M = K x M x N with periodic boundary conditions. Label sites of £ 
by indices {m,j, k}, with 1 < m < M, 1 < j < N and 1 < k < K. Label 
the state of the horizonal (vertical) edge incident at the site {m,j, k} in the 
direction of, say, decreasing {m,j} by a m jk ((3 m jk)- Denote the vertex weight 
at site {m,j, k} by W m jk and the interacting Boltzmann factor between the 
{m, j} sites of two adjacent layers k and k + 1 by B m j k . It can then be 
verified that the interaction given in Table 1 can be written in the form of 

B mjk = e W [Ka^-a 3+l ^)\ (16) 

where for convenience we have suppressed the subscripts m and k and intro- 
duced the notation 

Pm+l,j,k y fij, Pm,j,k+1 ~ > Pji (17) 

and similarly for the a's. Note that as a consequence of the ice rule the 
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quantity 



i JV i N 

i=i i=i 



is conserved within each layer and is independent of m. 

The problem at hand now is the evaluation of the partition function 

K M N , v 

zmnk = e e n n u( B mjkW mjh ) (19) 

"mjfe /3 mjfe fe=l m=l j=l V 

where the summations are taken over all edge states a m j k and [3 m jk, and the 
per-site "free energy" 

/= lim (MNK)" 1 In Z M nk- (20) 

M,N,K^oo 

3.2 The transfer matrix 

The partition function (|T^) can be evaluated by applying a transfer matrix 
in the vertical direction. In a horizontal cross section of C there are NK 
vertical edges. Let {(3 m } = {(3 mjk \l < j < N, 1 < k < K}, 1 < m < M 
denote the states of these NK vertical edges, and define a 2 NK x 2 NK matrix 
T with elements 

K N 

T({(3 m },{f3 m+1 }) = E II T[(B mj kW mjk ). (21) 

a mjk k=l j=l 

Then one has 

M 

Zmnk = E II T ({P m }, {Pm + i}) ~ A^ ax , (22) 

where A max is the largest eigenvalue of T. The validity of the ice rule within 
each layer now permits one to use a global Bethe ansatz in the analysis. 



The interlayer interaction (ITO) leads to a considerable simplification of 



the transfer matrix. It can be shown 11 that we have 



M N M 

II II B mjk = II ex P Nh(a mXk y k+ i - a mXk+ iy k ) . (23) 

m=l j=l m=l 
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As a result, only the conserved quantities y k and the state a m ^ k appear in 
the final product. This leads us to recast the partition function fll9|) in the 
form of 

M 

^MNK = 

E UT CS ({(3 m },{Pm + l}) 

= Tr(T eff ) M (24) 

with 

K , N , 

T eS ({P m }, {Pm+i}) = E II (e Nha -^ y ^-y^ I] Wmjk). (25) 

Umjk fc=l ^ j = l ' 

The problem is now reduced to one of finding the largest eigenvalue of the 
matrix T cfr with matrix elements as shown. The task is now considerably 
simpler since one needs only to keep track of two-dimensional systems. 



3.3 The free energy 



The eigenvalues of the effective matrix are obtained by applying a global 
Bethe ansatz consisting of the usual Bethe ansatz for each layer. The algebra 
is straightforward and the result is 

K 



hiNK ~ max TT {A L (n k )Y 

{n k } k=1 



with 



A L (n k ) = eMn^-n^^N-n, J| 



(k) 



UJ4 



(26) 



(27) 



being the eigenvalue for a m 1 k 

(k) 



—1 and, for each 1 < k < K, the n k complex 



numbers Za are 



Ik) 



e J e 



1 3 

\n h +l 



l,2,---,n k 



(28) 



where e j are n k distinct iVth roots of (— l) nfc+ . Using (PI) and (|20D , this 
leads to the per-site free energy 

1 i, 1 r< l -vk)n 

f = In u + lim max — 



In u + lim max , , 

K-+00 -l<y k <l K ^ 27T J-TT(l-y k )/2 



X 



\u u J 



(29) 
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3.4 The phase diagram 



Analyses of the free energy ( |29| ) lead to the following: ||11|| For large u, v or w, 
the system is frozen with complete ordering of u, v, or w dimers in all layers, 
and hence the free energies 

fu = Inu, U phase 

fv = hiv, V phase 

fw = lniy, W phase. (30) 

These are referred to as the U, V, and W phases, respectively. For large h, 
it is seen from Table 1 that the energetically preferred state is the one in 
which each layer is occupied by one kind of dimers, u, v, or w, and that the 
layers are ordered in the sequence of {u, w, v , u, w, v ■ ■ •}. This ordered phase 
is referred to as the H phase with the free energy 

fu = - \n{uvwe ih ), H phase. (31) 



More generally one finds that the extremum value in (g9|) always repeats 
in multiples of 3, namely, with 

Vk+3 = Vk- 

The following extremum sets of {yk} are found: 

1. {yi, jj2, 2/3} = {1, 1, 1}: In this case we have all y& = 1, and hence from 

(H 

f = fu, U phase. (32) 

2. {yi, 1/2,1/3} = { — — 1, — 1} : In this case we have all y& = —1, and hence 
from ((291) 



1 f n 

f = In u H — / In 

7T Jo 



- + -e td 
u u 



de 



= fw, w > v W phase 

= fy, v > w. V phase (33) 

3. {yi, y 2 , 2/3} = {1,-1, —1}: Substituting this sequence of yu values into ( p9|) 
and making use of the integration formula (0) in the resulting expression, one 
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obtains 



1 f w 

f = lntt + — / In 

D7T J -TV 


u u 


1 r 
d6 + — In 

07T J -TV 


E + V 2 V* 
u u 




—Ah Ah 

ve < ve < w 






w < ve 


4h < ve 4h 




= Sh- 


—Ah Ah 

ve <w< ve . 





dO 



(34) 

Now the free energies in the second and third lines can be discarded since 
they are always smaller than the largest of {fu, fv, fw}- Thus, this set of 
{i/k} leads to a frozen ordering for sufficiently large h as indicated in the last 
line. 

4. {ui, y 2 , 2/3} = {y, y, y}: In this case all y^ = y, y maximizes the free energy 
(|29|). Then, substituting y k = y into p9|) and carrying out a straightforward 
differentiation with respect to y, one obtains 



where the extremum yo is given by 



Y phase 



2(l-Z/0) 



r U 2 



COS 



ur 



,2n 



2wv 



(35) 



(36) 



This is a disorder phase which we refer to as the Y phase. 

5. {yi, 2/2, 1/3} = {yi, —1,-1}: This is the H phase with the u layers replaced 
by layers with y k = y 1 , so that the layer ordering is {yi, w, v, yi, w, v • • •}. 
This is a partially ordered phase which we refer to as the I u phase. 

6- {2/1,^/2,2/3} = {1,2/2, _ 1} : This is the H phase with the w layers replaced 
by layers with y k+ i = y 2 , so that the layer ordering is {u,y 2 ,v,u,y 2 ,v ■ ■ ■}. 
We refer to this as the I w phase. 

7. {2/1,2/2,2/3} — {1, — 1, 2/3}: This is the H phase with the v layers replaced 
by layers with 2/^+2 = 2/3, so that the layer ordering is {u, w, 2/3, u, w, 2/3 • • ■}. 
We refer to this as the /„ phase. 

Since the phase diagram must reflect the {it, v, w} symmetry of the inter- 
layer interaction given in Table 1, it is convenient to introduce coordinates 

X = hx(v/w) Y = ln(i;ii;/n 2 ) (37) 
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such that any interchange of the three variables u, v, and w corresponds to a 
120° rotation in the {X, Y} plane. The phase boundaries are then determined 
by equating the free energies of adjacent phases. The results are collected in 
Fig. 3. One finds the following: 

• h < h : For small h the phase diagram is the same as that of the 
h — 0, the diagram shown in Fig. 3(a). The phase boundaries between the 
{U, V, W} phases and the Y phase are 

u — \v ± w\. (38) 

• ho < h < hi. As h increases from zero, our numerical analyses indicate 
that the H phase appears when h reaches the value 

3 /-27T/3 

h = — ln(2 + 2 cos 6)de = 0.2422995.... (39) 
8n Jo 

The resulting phase diagram is shown in Fig. 3(b). The phase boundary 
between the H and Y phases is 

Un(uvwe 4h ) = f Y (y ). (40) 

• hi < h < h,2- As h increases from ho, the I u ,Iv, Iw phases appears when 
h reaches a certain value h\. The resulting phase diagram is shown in Fig. 
3(c). The numerical value of hi is given by 

h x = - ln(—) = 0.2552479 
4 V u ) 

where v ju is the solution of the two equations 

V 1 1 /-t(1-J/o)/2 

(l+y )ln- + -(l + 3y )m2 = -/ ln(l + cos0)d0 (41) 

U D 7T J0 

= cos "(^- 1 )- (42) 

• hi < h < h%: As h increases from hi, it was found that the regimes I u , 
I v , and I w extends to infinite when h exceeds the value 

h 2 = (ln3)/4 = 0.2746530- ••. 
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The phase diagram is shown in Fig. 3(d). 

• /13 < h < h±. As h increases from hi, it was found that the boundary 
of the H phase bulges toward the U, V, W phases along the 30°, 150°, 270° 
lines, touching the U, V, W boundaries in these directions when h equals to 

h 3 = (ln2)/2 = 0.3465735- •-. 

For h > h 3 , the H phase borders directly with the U, V, W phases with 
respective boundaries 

u 2 = vwe 4h , v 2 = wue* h , w 2 = uve 4h . (43) 

The size of these borders grows while the Y phase shrinks as h increases. 
The phase diagram in this regime is shown in Fig. 3(e). 

• h > hi'. As h increases further from h^, it was found that the Y phase 
disappears completely when h exceeds the value 

h A = - ln(l + 7 2 ) = 0.3816955 • • • , 
8 

where 7 is the solution of the equation 

7 — tan^ 1 7 = 7r. (44) 



3.5 The critical behavior 

The critical behavior near all phase boundaries can be analyzed by expand- 
ing the free energy. We find all transitions to be of first order, except 
those between the {U,V,W} and Y phases, and between the {I U ,I V ,I W } 
and H phases, which are shown below to be the same as that in the 5- vertex 
model, [|n| |TB| namely, a continuous transition with a square-root singular- 
ity in the specific heat. This transition, first reported by this author f|T~6| 



in 1967, is now more widely known as the Pokrovsky-Talapov type phase 



transition. 17 



To see the occurrence of the Pokrovsky-Talapov type transitions, we 
note that the ordered U, V, W and H phases (with y k = ±1) have constant 
free energies and hence zero first derivatives. The transition between the 
U, V, W phases and the Y phase is the same as in two dimensions, and is 
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continuous. ||15|| This fact can be seen by expanding the free energy near %/q 
as 



fy{y) = fr(yo) + {y- yo)f Y (yo) + h{y - yo) 2 f Y (yo) 
+^(y-yoff Y '(yo) + ---. 



(45) 



Using the expression of f Y (y) given by (|35D, one sees that, indeed, the first 
derivative f Y (yo), yo — ±1, vanishes identically on the boundary (|38"D which 
is precisely f Y (yo) = 0. Furthermore, it is also seen that f Y (yo) ~ sin[7r(l — 
2/o)/2] = 0. Therefore, the extremum of f Y (y) given by ([45]) occurs at y = 
2/cxtrm determined from 



|/extrm — I/O — =ti 



2f Y (y ) 

\ -mvo) 



(46) 



where t = \T — T c \, T c being the critical temperature. Substituting this j/ ex trm 
into one obtains 



/y(2/extrm) 



friyo) ± g/y(2/o 



2Myo) 



f Y (y Q ) + c(u,v,w,h)t 3/2 . 



(47) 



This leads to a square- root singularity in the specific heat. 

Applying the same analysis to the H and I u phases, one concludes 



extrm ) = /(2/10) + c i( u j v > w > h)t 3 ^ 2 . (48) 
This gives rise again to a square- root singularity in the specific heat. 



4 Summary 

We have described the exact solutions of two lattice models in three di- 
mensions, an 0{— 1) model solved using the combinatorial method of Pfaf- 
fians, and a layered dimer model solved using an algebraic transfer matrix. 
However, both models stop short of being genuinely three-dimensional. The 
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0(—l) model suffers two defects. First, it describes line configurations run- 
ning only in preferred directions. Secondly, the Boltzmann weights can be 
negative. The solution of a model without these restrictions would help 
us to solve the Ising model in three dimensions. The layered dimer model, 
while having strictly positive weights, describes dimer configurations in which 
dimers are confined in planes. As a consequence, the critical behavior is es- 
sentially two-dimensional. 
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Figure captions 

Fig. 1. (a) A lattice point of the simple cubic lattice C. (b) The split of the 
lattice point in (a) into two points. 

Fig. 2. A stack of K layers of honeycomb dimmer lattices. 
Fig. 3. The phase diagram. 
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Fig. 2 
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Fig. 3 



